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Abstract 

Solid  oxide  fuel  cell  (SOFC)  is  a  complicated  system  with  heat  and  mass  transfer  as  well  as  electrochemical  reactions.  The  real-time  dynamic 
simulation  of  SOFC  is  still  a  challenge  up  to  now.  This  paper  develops  a  one-dimensional  mathematical  model  for  direct  internal  reforming  solid 
oxide  fuel  cell  (DIR-SOFC).  The  volume-resistance  (V-R)  characteristic  modeling  technique  is  introduced  into  the  modeling  of  the  SOFC  system. 
Based  on  the  V-R  modeling  technique  and  the  modular  modeling  idea,  ordinary  differential  equations  meeting  the  quick  simulation  are  obtained 
from  partial  differential  equations.  This  model  takes  into  account  the  variation  of  local  gas  properties.  It  can  not  only  reflect  the  distributed  parameter 
characteristics  of  SOFC,  but  also  meet  the  requirement  of  the  real-time  dynamic  simulation.  The  results  indicate  that  the  V-R  characteristic  modeling 
technique  is  valuable  and  viable  in  the  SOFC  system,  and  the  model  can  be  used  in  the  quick  dynamic  and  real-time  simulation. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  energy  conversion  devices 
that  produce  electricity  and  heat  directly  from  a  gaseous  or  gasi¬ 
fied  fuel  by  electrochemical  reaction  with  an  oxidant  [1].  SOFC 
operates  at  high  temperature  and  atmospheric  or  higher  pressure, 
and  can  use  hydrogen,  carbon  monoxide,  and  hydrocarbons  as 
fuel,  and  air  or  pure  oxygen  as  oxidant.  As  a  matter  of  fact,  both 
simple-cycle  and  hybrid  SOFC  systems  have  been  demonstrated 
among  the  highest  efficiencies  of  any  power  generation  system, 
combined  with  minimal  air  pollutant  emissions  and  low  green¬ 
house  gas  emissions.  These  capabilities  have  made  SOFC  an 
attractive  emerging  technology  for  stationary  power  generation, 
especially  for  the  distributed  generation  [2] . 

The  conventional  modeling  technique  for  a  distributed 
parameter  SOFC  is  to  establish  a  set  of  partial  differential  equa¬ 
tions  based  on  the  conservation  equations  for  species  mass,  total 
mass,  momentum  and  energy.  The  solution  to  these  equation  sets 
is  very  difficult  and  very  time  consuming.  In  addition,  sometimes 
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the  exact  numerical  solution  may  be  impossible.  However,  the 
SOFC  dynamic  and  real-time  simulation  is  of  significance  in  the 
verification  of  controller  or  online  monitoring  system.  In  order 
to  overcome  this  problem  and  also  provide  reference  for  opti¬ 
mization  design  of  the  SOFC,  we  need  a  real-time  model  which 
is  also  a  time  saving  model. 

The  main  purposes  of  this  paper  are  as  following:  first, 
considering  the  variation  of  local  flow  properties,  such  as  pres¬ 
sure,  mass  flow  rate,  density,  specific  heat  capacity,  thermal 
conductivity  and  dynamic  viscosity;  secondly,  introducing  the 
volume-resistance  characteristic  modeling  technique  into  the 
SOFC  system  to  establish  a  modular  model  satisfying  the  quick 
dynamic  simulation;  finally,  developing  a  dynamic  model  of  a 
planar  DIR-SOFC  cell. 

2.  SOFC  mathematical  model 

There  have  been  several  publications  focusing  on  modeling 
the  performance  of  SOFCs.  Such  models  can  be  for  differ¬ 
ent  geometries,  PEN  (positive-electrode/electrolyte/negative- 
electrode)  structures,  flow  configurations,  and  operating 
temperature  ranges,  etc.  [1,3-20].  Most  publications  have  been 
addressed  based  on  CFD  with  detailed  parameters  or  steady  and 
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Nomenclature 

Aj  cross-section  area  of  the  channel  (m2) 

Cp  specific  heat  capacity  (kJ  kg- 1  K- 1 ) 
d\y  hydraulic  diameter  (m) 

e  specific  internal  energy  (kJ  kg- 1 ) 

/  resistance  coefficient 

F  Faraday  constant  (C)  or  molar  flow  rate  (mol  s- 1 ) 

G  mass  flow  rate  (kg  s-1) 
hj  height  of  the  channel  (m) 

H  specific  enthalpy  (kJ  kg  - 1 ) 

7,  7  local  and  average  current  density  (A  m-2) 

£/,pen>  kp intc  air  and  fuel  channels  heat  transfer  coeffi¬ 
cient  (kJ  m-1  s-1  K-1) 

L  width  of  the  activation  area  of  the  fuel  cell  (m) 

LHV  lower  heating  value  (kJ  mol-1) 

P  pressure  (Pa)  or  power  density  (W  m-2) 

p  partial  pressure  (Pa) 

R  gas  constant  (J  mol- 1  K- 1 ) 

Uj  velocity  in  gas  channels  (ms-1) 

U  voltage  (V)  or  utilization 

W  width  of  the  activation  area  of  the  fuel  cell  (m) 
yi  molar  fraction  of  component  i 

Greek  letters 

AH  enthalpy  change  of  reaction  (kJ  mol-1) 

y)  fuel  electrical  efficiency 

ka,  Xf  thermal  conductivity  of  the  air  and  fuel  steam 

(kJ  m-1  s-1  K-1) 

^pen>  ^intc  thermal  conductivity  of  PEN  and  intercon¬ 
nector  (kJ  m-1  s-1  K-1) 
fi  dynamic  viscosity 

Vifk  stoichiometric  coefficient  of  component  i  in  reac¬ 
tion  k 

p  density  of  gas  steams  and  solid  parts  (kg  m-3) 

a  Stefan-Boltzmann  constant  5.676E— 8 

(W  m-2  K-4) 

^anode?  ^cathode  electronic  conductivity  of  anode  and  cath¬ 
ode  (£2-1  m-1) 

g electrolyte  ionic  conductivity  of  electrolyte  (Q~l  m-1) 
r  thickness  of  solid  structures  (m) 

G  emissivity 

Subscripts 
a  air  channel 

anode  anode  polar 

cathode  cathode  polar 

electrolyte  electrolyte 
f  fuel  channel 

INTC  interconnector 

ocp  open  circuit  potential 

PEN  PEN  structure 

Superscripts 

in  fuel  and  air  inlet 


Interconmctor 

CHA,H2tCO, 

[l)CHA+H2Oc*CO+3H2 

C02,H20 

[1I)C0+H20c*C02  +  H2 

■v  h2  h& 

(HI)H,  +0:~  H:0+ 2«-\  / 

Anode 

to2- 

Electrolyte 

(jp)o.5a+2e"-»a’" 

Cathode 

(X.N, 

j 

Interconnect 

CHi,HlrCO, 
CO,,  H.0 


o>.n2 


Fig.  1.  Schematic  principle  view  of  a  co-flow  planar  SOFC  cell. 


dynamic  model  with  many  assumptions  and  simplifications,  for 
example  the  flow  properties  and  gas  velocity  are  taken  as  con¬ 
stant  throughout  the  system  based  on  the  inlet  conditions.  In 
recent  years,  several  papers  (Iora  et  al.  [3],  Janardhanan  et  al. 
[19],  Liu  et  al.  [20])  have  considered  the  variation  of  local  gas 
properties.  But  owing  to  the  complicated  numerical  iteration, 
these  models  could  not  meet  the  requirement  of  the  real-time 
dynamic  simulation. 

The  dynamic  SOFC  model  developed  here,  which  is  one¬ 
dimensional  planar  DIR-SOFC  model,  considers  the  variation 
of  local  flow  properties,  and  allows  for  co-flow  and  counter¬ 
flow  operations.  The  approach  to  avoid  the  iteration  of 
the  coupling  of  pressure  and  mass  flow  rate  is  based  on 
the  volume-resistance  characteristic  modeling  technique,  the 
distributed-lumped  parameter  method  and  the  modular  mod¬ 
eling  idea.  The  established  model  can  not  only  reflect  the 
characteristics  of  distributed  parameter,  but  also  satisfy  the 
requirement  of  the  transient  simulation.  Fig.  1  presents  the 
schematic  principle  view  of  a  co-flow  planar  DIR-SOFC  cell. 

To  simplify  the  modeling,  some  reasonable  assumptions  are 
summarized  as  follows: 

(1)  Only  H2  is  electrochemically  oxidized. 

(2)  Constant  Nusselt  number. 

(3)  Adiabatic  boundaries  for  the  cell. 

2.7.  Species  mass  balances  model 

The  chemical  species  considered  in  the  fuel  channel  are  CH4, 
H2,  CO,  CO2  and  H2O,  while  in  the  air  channel  the  chemical 
species  are  O2  and  N2  (as  the  oxidant  is  considered  to  be  air 
and  not  pure  oxygen).  In  Table  1  the  reactions  and  the  corre¬ 
sponding  reaction  rates  in  the  SOFC  system  are  presented.  It  is 
assumed  all  CO  is  converted  through  the  water  gas  shift  reac¬ 
tion  (II),  which  is  at  equilibrium;  all  CH4  in  the  fuel  is  only 
concerned  with  the  steam  reforming  reaction  (I).  The  molar  flux 
in  both  gas  channels  is  considered  convective  in  the  flow  direc¬ 
tion.  It  is  reasonable  to  neglect  the  axial  dispersion  effects  [1]. 
The  detailed  gas  properties  (density,  specific  heat  capacity,  ther¬ 
mal  conductivity,  dynamic  viscosity,  etc.)  are  mainly  affected 
by  the  local  temperature  and  the  components,  their  detailed  cal¬ 
culations  can  be  found  in  the  references  [22,23].  For  instance, 
the  dynamic  viscosity  of  multi-component  gas  mixtures  is  based 
on  the  Reichenberg’s  expression,  and  the  thermal  conductivity 
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Table  1 

Reactions  and  reaction  rates  in  the  SOFC  system 


Steam  reforming  [1,4,5] 

(I)CH4  +  H20  O  CO  +  3H2 

P(I)  =  &act/?f,CH4  exp 

f  ^act  \ 

\~RTf) 

(1) 

Water  gas  shift  [1,4] 

(II)C0  +  H20<S>C02  +  H2 

P(Il)  =  &WGSRPf,CO  ^ 

1  Pf,co2Pf,n2  \ 

^shift  Pi,  CO  Pi,  H2  0  J 

(2) 

Oxidation 

Reduction 

(III)H2  +  O2-  -»•  H20  +  2e-1 
(rV)0.5O2  +  2e_1  ^02“ 

I 

R(m,w)  =  Yp 

(3) 

Table  2 

Species  mass  balances  equations  of  SOFC  model 

Fuel  channel :  — ^  = 
dt 

_dCfi!Uf  +  I H.kRkf  i  e  {Cftt,  H2,  CO,  CO2,  H20)  (4) 

dx  7-^  hf 

k£  {(I), (II), (III)} 

3  Caz-  3  Ca;wa  1 

Air  channel  :  -—f-  = - -1 - h  (iv)  7- ,  i  e  {02,  N2}  (5) 

cue  /z  ^ 


of  multi-component  gas  mixtures  is  based  on  the  Wassiljewa’s 
expression  and  the  Mason  and  Saxena  modification  (Poling  et 
al.  [23]). 

The  kinetics  of  the  reforming  reaction  (I),  the  water  gas 
shift  reaction  (II)  and  the  electrochemical  reaction  (III),  which 
are  listed  in  Table  1.  The  first  order  kinetic  expression 
derived  by  Achenbach  and  Riensche  [5]  has  been  adopted 
in  the  reforming  reaction  (I),  with  an  activation  energy 
of  £’act  =  82x  103Jmol-1  and  a  pre-exponential  constant  of 
&act  =  0.04274  mol  s-1  m-2  Pa-1 .  Such  formation  is  considered 
typically  for  the  DIR-SOFC  performance  [1].  The  water  gas 
shift  reaction  (II)  is  considered  a  very  fast  one  and  is  assumed  to 
be  at  equilibrium  in  the  fuel  channel  [6] .  The  equilibrium  lim¬ 
ited  shift  reaction  rate  expression  is  derived  by  Aguiar  et  al.  [1], 
where  &shift  =  exp(4276/7f  —  3.961)  [7].  The  kinetic  of  the  elec¬ 
trochemical  reaction  is  only  related  to  the  local  electric  current 
density. 

The  detailed  species  mass  balances  equations  listed  in  Table  2 
are  referred  to  the  literature  [3]. 

2.2.  Mass  balances  model 


the  same  to  the  electrochemical  reaction  rate,  and  Mo2  is  the 
oxygen  molecular  weight. 

2.3.  Energy  balances  model 


In  a  SOFC  system,  current  and  temperature  distributions  are 
strongly  coupled.  Therefore,  exact  energy  balances  model  and 
accurate  temperature  profiles  in  the  fuel  cell  are  very  impor¬ 
tant  and  essential.  The  model  in  Table  4  considers  that  the 
thermal  flux  in  both  the  PEN  structure  and  interconnector  is 
conductive,  which  is  modeled  by  Fourier’s  law  of  heat  con¬ 
duction.  Meanwhile,  the  thermal  flux  in  the  gas  channels  is 
convective  from  the  gas  channels  to  the  solid  parts  [1].  The 
heat  transfer  coefficients  between  the  gas  channels  and  the 
solid  parts  were  calculated  using  a  constant  Nusselt  number 
(here,  it  is  taken  to  be  3.09;  Aguiar  et  al.  [1]).  In  general, 
an  assumption  is  given  for  the  laminar  flow  conditions  that 
the  Nusselt  number  is  considered  independent  of  the  Reynolds 
number  [24] .  Due  to  the  high  temperatures  involved,  radiation 
energy  between  the  PEN  structure  and  interconnector  is  also 
included. 

Here,  the  following  enthalpy  change  in  every  reaction  is 
derived  according  to  the  method  depicted  detailed  in  the  ref¬ 
erence  [25],  where  A H°  is  the  enthalpy  change  at  the  standard 
state: 

A  Hi  =  AH°i  -  16373.61  +  R  ( 7.9517f  -  4.354e  -  3rf2 


+  0.7213e  —  6Tf 


0.097e5\ 
Tf  ) 


(14) 


For  the  mass  balance,  oxygen  ions  produced  in  the  cathode 
transfer  to  the  PEN  structure.  Hydrogen  in  the  fuel  channel  will 
penetrate  into  the  PEN  structure  as  well,  where  the  electrochem¬ 
ical  reaction  occurs.  The  water  steam  produced  from  the  reaction 
transfers  to  the  fuel  channel.  Therefore  the  increased  mass  in  the 
fuel  channel  is  only  owing  to  the  mass  transfer  of  the  oxygen 
ions.  For  the  whole  system  the  mass  is  conserved.  The  detailed 
Eqs.  (6)  and  (7)  are  listed  in  Table  3,  where  vo2  is  stoichiometric 
coefficient  of  oxygen,  Rq2  the  oxygen  reaction  rate,  in  fact  it  is 


AHn  =  AH°u  -  7756.56 


+  R 


^1.867f-0.27e-37’f2  + 


1.164e5\ 
Tf  ) 


(15) 


0.1515e5\ 


A  Hm  =  AH°m  +  4097.22 

+  R  ( —  1.59857pEN+0*3875£  —  3TpEN - 1 

Tpen  ) 

(16) 


Table  3 

Mass  balances  equations  of  SOFC  model 

dpf  dpfUf  1 

Fuel  channel  :  —  = - vo?  Ro0  Mo?  —  (6) 

dt  dx  hf 

dp  PL  dp  pi  Upl  1 

Air  channel  :  —  = - - - h  vo2Ro2Mo2  —  (7) 

dt  dx  h  a 


2.4.  Momentum  balances  model 

To  calculate  the  pressure  profiles  in  the  fuel  and  air  channels,  a 
momentum  balance  for  each  gas  channel  is  required.  The  energy 
conservational  equations  of  Iora  et  al.  [3]  revised  a  little  are 
adopted,  which  are  presented  in  Table  5. 


Table  4 

Energy  balances  equations  of  SOFC  model 
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Table  5 

Momentum  balances  equations  of  SOFC  model 


Fuel  channel  : 


Air  channel  : 


dpfUf 

dpfUfUf 

dE 

2  f(W,  hf)  •  utiLf 

'  dt  = 

dx 

dx 

d2 

dpaUa 

dpaUaUa 

dPa 

2 /(W,  ha )  •  wa/za 

dt 

dx 

dx 

d2 

(17) 

(18) 


2.5.  Electrochemical  model 

The  electrochemical  model  here  considers  the  distribution 
parameter.  The  open  circuit  potential  Uoc p  is  determined  by  the 
gas  composition  in  both  gas  channels  and  local  temperature  at  the 
PEN  structure.  It  is  depicted  by  the  Nernst  equation  as  follows: 


Uocp  =  <  + 


h2 


RTpEN  ^  (  Pf,H2(Pa,02)' 

2  F 


,0.5 


In 


Pf,H2o 


(19) 


(20) 


R7pen  i  )_L) 

4  F  \P^J 

where 

u^2  =  1.2723  -  2.7645e  -  4  x  7pEN 

C/h2  is  the  ideal  voltage  for  hydrogen  oxidization  [2,8].  Pst(i 
is  taken  as  the  standard  pressure. 

Only  activation  and  ohmic  losses  are  taken  into  account,  and 
concentration  loss  is  neglected  [1].  Ohmic  resistances  are  caused 
by  resistance  to  conduction  of  ions  and  electrons  and  contact 
resistance  between  cell  components  [1,9,10,11]: 


D  7 mode  Te]ectrolyte  Tcathode 

Tvohm  —  I  I 

^anode  ^electrolyte  ^cathode 


(21) 


-2 


Overpotentials  at  the  cathode  and  at  the  anode  are 
assumed  independent  of  the  local  current  density.  The 
model  derived  by  Achenbach  [6]  is  adopted,  where  activa¬ 
tion  energies  Eanode  =  1 1 0  kJ  mol- 1 ,  ^cathode  =  1 60  kJ  mol- 1 
and  corresponding  coefficients  &anode  =  2.13E+8  Am 

^cathode  : 

1 

^anode 

1 


1.49E+10  Am-2,  m  =  0.25: 


:  ^anode 


RTf  V  Pf  ) 


(Fanode/^7f) 


(22) 


4  F  f  pp2 

kathode  '  Calh0de  RTa  V  Ai 


(^cathode/  PTf) 


(23) 


As  both  electrodes  are  normally  good  conductors,  a  con¬ 
stant  cell  operating  voltage  Uop  throughout  the  cell  is  normally 
considered  [1,14]. 

Some  fuel  cell  performance  parameters  are  defined  as  fol¬ 
lows: 


Current  density  :  I  = 


Fuel  utilization  :  Uf  - 


U ocp  U op 


Air  utilization  :  U*  = 


^ohm  T  ^anode  H-  ^cathode 
IWL 

:2F(4>gl4  +  3>S2+3&  )lf 

IWL 


4  FySF? 


(24) 


(25) 


(26) 
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Power  density  :  P  =  IU 


(27) 


Fuel  electrical  efficiency  :  r\ 

PLW 

~  ^h4LHV^H4  +  yj?2LHVg2  + 


(28) 


3.  SOFC  V-R  mathematical  model 


3.1.  V-R  characteristic  modeling  technique 


It  is  difficult  to  solve  partial  differential  equations.  More 
importantly,  the  coupling  of  pressure  and  flow  rate  will  cause 
the  iteration  solution  to  the  equations.  Here  the  volume  and 
resistance  (V-R)  characteristic  modeling  technique  [21]  is  intro¬ 
duced  into  the  complicated  SOFC  system  modeling  to  establish 
ordinary  differential  equations  model.  Some  state  variables  were 
introduced  into  the  model  to  decouple  the  pressure  and  flow  rate. 

In  the  fluid  system,  each  component  can  be  treated  as 
one  of  three  types:  volume  module,  resistance  module  and 
volume-resistance  module.  The  volume  module  is  the  compo¬ 
nent  which  can  neglect  the  pressure  loss,  and  accumulate  fluid 
based  on  the  difference  of  flow  rates  between  the  inlet  and  the 
outlet.  The  resistance  module  is  the  component  where  the  fluid 
has  an  obvious  pressure  loss,  i.e.  the  pressure  of  the  inlet  and 
outlet  is  different,  which  determines  the  flow  rate  of  the  compo¬ 
nent.  To  some  extent,  many  fluid  components  have  the  qualities 
of  both  the  volume  module  and  the  resistance  module.  The  SOFC 
system  is  just  this  kind  of  component. 

The  V-R  characteristic  modeling  technique  is  based  on  the 
lumped  parameter  modeling  idea. 

In  general,  the  volume  module  can  be  described  as 


d  P 
d  t 


RT 

"chT 


(Ci  -  G2) 


(29) 


On  the  other  hand,  the  resistance  module  can  be  described  as 


dG 
d  t 


A 

-  —  (Pi-P\)- 

dx 


2/G2 

Apdh 


(30) 


In  order  to  depict  the  distributed  parameter  SOFC  cell,  several 
control  units  must  be  divided  in  the  SOFC,  each  unit  will  be 
treated  as  a  lumped  parameter  module.  In  each  unit,  the  module 
will  be  separated  into  two  sub-modules.  One  is  a  volume  sub- 
module  and  the  other  is  a  resistance  sub-module. 

Fig.  2  illustrates  how  the  volume-resistance  characteristic 
modeling  technique  is  used  in  the  flow  grid.  In  the  conventional 
modeling  method,  the  pressure  and  flow  rate  for  the  control  unit 
should  be  solved  using  iteration  method.  While  in  the  V-R  mod¬ 
eling  method,  the  control  unit  will  be  divided  into  two  modules. 
In  the  volume  module,  the  pressure  is  considered  constant,  and 
the  mass  flow  rate  is  the  state  variable.  In  the  resistance  module, 
the  mass  flow  rate  is  constant,  and  the  pressure  is  the  state  vari¬ 
able.  Thus,  it  is  no  necessary  to  assume  pressure  to  solve  flow 
rate  through  iteration  algorithm,  the  pressure  and  flow  rate  will 
be  decoupled  through  two  state  variables  in  individual  modules. 

The  key  point  of  the  VR  model  is  non-iterative  algorithm  in 
the  calculation.  Therefore,  it  will  save  calculation  cost,  and  the 
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Fig.  2.  Diagram  of  grid  mark  for  mass  flow  rate  and  pressure  in  one  section. 


most  important  point  is  that  calculation  time  can  be  determined. 
While  for  the  conventional  iteration  algorithm,  the  convergence 
time  cannot  be  determined  owing  to  the  uncertainty  of  iteration 
times.  The  non-iterative  algorithm  will  guarantee  VR  model  to 
be  the  real-time  model,  which  can  be  used  in  the  online  simula¬ 
tion. 

3.2.  SOFC  V-R  characteristic  model 

Based  on  the  above  idea,  one  set  of  equations  for  the  fluid  in 
the  control  section  can  be  obtained  in  Table  6. 

Here  for  the  gas  channels,  using  the  mass  balance  theorem, 
the  following  Eqs.  (32)  and  (36)  can  be  obtained  from  the  volume 
sub-module;  using  the  momentum  balance  theorem,  the  Eqs. 
(33)  and  (37)  from  the  resistance  sub-module  can  be  obtained. 
While  dealing  with  the  species  mass  balance  and  energy  balance, 
the  two  sub-modules  are  merged  into  one  control  section. 

For  the  solid  parts,  the  energy  balance  equations  of  both  the 
PEN  structure  (12)  and  the  interconnector  (13)  are  considered 
in  one  control  section. 

4.  Solution  strategy 

4.1.  Boundary  and  initial  conditions 

The  geometry  parameters  should  be  given  based  on  the 
defined  fuel  cell  stack.  All  the  static  variables  should  be  given  an 
initial  value  to  start  the  simulation,  such  as  molar  components, 
molar  flow  rates,  temperatures  of  the  fuel  channel  and  the  air 
channel  as  well  as  the  PEN  structure  and  interconnector,  etc. 

At  the  same  time,  boundary  condition,  i.e.  the  thermal  flux 
density  of  both  the  inlet  and  the  outlet  for  the  solid  parts  are 
considered  as  zero.  In  addition,  the  pressure  of  both  channel 
outlets  is  given  as  certain  value.  As  a  result,  compared  with  the 
middle  modules,  the  difference  algorithms  of  both  the  front  and 
the  end  modules  are  different.  So  the  modules  at  the  front  and 
the  end  are  treated  independently. 

4.2.  Simulation  method 

After  the  conversion  from  the  partial  differential  equa¬ 
tions  to  the  ordinary  differential  equations  based  on  the 
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SOFC  volume-resistance  characteristic  model 
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State  variables  (Num:15) 
Molar  concentration  ;  C 
Fuel  channel: 

ch4,h2,co,co2  h2o 

Air  channel: 

o2,n2 

Mass  flow  rate  :  G 
G-air,  G-fuel 

Pressure :  P 
P-air,  P-fuel 

Temperature  :  T 
Air  channel.  Fuel  channel, 
Interconnect,  PEN 

Fig.  3.  Flowing  interface  of  the  control  section  in  SOFC  cell. 


distributed-lumped  parameter  method  and  the  volume- 
resistance  characteristic  modeling  technique,  the  single  module 
of  every  unit  is  given  in  Fig.  3. 

The  system  model,  which  is  based  on  the  above  mathematical 
model,  was  developed  on  the  EASY5  simulation  platform  using 
the  FORTRAN  language  compiler. 

5.  Result  and  discussion 

In  this  section,  a  specific  co-flow  fuel  cell  is  simulated.  Sec¬ 
tion  5.1  introduces  the  model  input  parameters  and  operating 
conditions;  Section  5.2  analyzes  the  steady-state  performance  of 
a  co-flow  SOFC  system;  Section  5.3  briefly  studies  a  transient 
process  with  a  step  change. 

5.1.  Input  parameters  and  operating  conditions 

The  geometry  input  parameters  and  some  other  physical  and 
electrical  properties  adopt  the  parameters  in  Ref.  [1].  The  main 
geometry  dimension  parameters  and  operating  conditions  are 
listed  in  Table  7. 


Table  7 

Key  geometry  parameters  and  operating  conditions 


Dimension  parameters  (unit:  m) 

Cell  length 

0.4 

Cell  width 

0.1 

Fuel  channel  height 

0.001 

Air  channel  height 

0.001 

Anode  thickness 

5E-4 

Cathode  thickness 

5E-5 

Electrolyte  thickness 

2E-5 

Interconnector  thickness 

5E-4 

Operating  conditions 

Inlet  fuel  and  air  temperature 

1173  K 

Outlet  fuel  and  air  pressure 

1.0  bar 

Inlet  fuel  molar  flow  rate 

0.0012  mol  s-1 

Inlet  air  molar  flow  rate 

0.012  mol  s-1 

Fuel  composition 

0.1710  CH4,  0.2626  H2,  0.0294  CO, 
0.0436  C02,  0.4934  H20 

Air  composition 

0.21  02,  0.79  N2 

Fuel  utilization 

0.85 

Fig.  4.  Fuel  and  air  channels,  interconnector  and  PEN  structure  temperature 
profiles  along  the  cell  length. 


5.2.  Steady -state  simulation  results 

The  steady- state  simulation  is  based  on  the  above  input 
parameters  and  operating  conditions. 

The  temperature  profiles  of  the  fuel  and  air  channels,  the  inter¬ 
connector  and  PEN  structure  along  the  cell  length  are  presented 
in  Fig.  4.  Owing  to  the  endothermic  reaction  of  the  reform¬ 
ing  of  methane  at  the  entrance,  the  temperature  decreases  by 
50  K.  Then  the  temperature  increases  along  the  fuel  and  air  flow 
directions  by  heat  accumulation,  which  is  released  by  the  elec¬ 
trochemical  and  water  gas  shift  reactions.  The  maximum  tem¬ 
perature  occurs  at  the  outlet.  In  general,  the  maximum  tolerable 
temperature  gradient  is  approximately  lOKcm-1  by  assuming 
a  maximum  safe  stress-induced  strain  of  0.1%  and  a  thermal 
expansion  coefficient  for  the  YSZ  material  [13].  Thus,  the  tem¬ 
perature  gradient  calculated  for  the  base  case  is  reasonable. 

In  this  type  of  fuel  cell  stack,  internal  reforming  reaction,  the 
water  gas  shift  reaction  and  the  electrochemical  reaction  will 
occur  simultaneously.  The  mole  fraction  profiles  along  the  cell 
length  of  all  the  species  in  the  fuel  channel  steam  are  illustrated 
in  Fig.  5.  Methane  will  be  reformed  rapidly  at  the  entrance  of 
the  fuel  cell.  Consequently,  the  concentration  of  the  hydrogen 
increases  fast.  Then  the  electrochemical  reaction  displays  the 
main  influence,  the  hydrogen  concentration  will  decrease  along 
the  cell,  the  carbon  dioxide  will  increase  correspondingly.  For 
the  rated  case,  at  the  exit  of  the  fuel  channel,  all  the  methane 
has  been  fully  consumed,  and  the  composition  is  7.97%  in  H2, 
2.93%  in  CO,  15.25%  in  C02,  73.84%  in  H20. 

Figs.  6  and  7  present  the  profiles  of  the  dimensionless  gas 
properties  (pressure,  velocity,  density,  specific  heat  capacity, 
thermal  conductivity  and  dynamic  viscosity)  along  the  cell 
length  in  the  air  and  fuel  channel,  respectively.  The  dimension¬ 
less  gas  properties  are  defined  to  be  the  ratio  of  local  values 
and  corresponding  inlet  gas  properties  parameters.  The  inlet  gas 
properties  can  be  calculated  by  the  inlet  conditions  using  the 
methods  depicted  in  Refs.  [22,23].  It  is  clearly  shown  that  the 
variation  of  the  gas  properties  is  lower  in  the  air  channel  than  in 
the  fuel  channel  [3]. 
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Gas  properties 
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The  pressure  drops  by  3.06%  and  0.12%  in  the  air  and  fuel 
channels,  respectively.  In  general,  the  flow  rate  and  temperature 
will  increase  along  the  cell  length  in  the  fuel  channel.  So  the 
pressure  should  increase  if  the  flow  resistance  is  not  considered. 


Fig.  8.  Mass  flow  rate  profiles  along  the  cell  length. 


Here,  owing  to  the  increasing  of  the  velocity,  the  pressure  loss 
will  increase  along  the  cell  length.  Therefore,  the  pressure  drop 
due  to  flow  resistance  will  be  dominated,  and  this  causes  the 
pressure  drop  in  the  fuel  channel. 

The  density  is  related  to  the  pressure  and  the  temperature 
based  on  the  gas  state  equation.  In  the  air  channel,  the  pressure 
has  little  change  at  the  entrance,  while  the  temperature  drops, 
so  the  density  in  the  air  channel  will  rise  at  the  entrance.  With 
the  increase  of  the  temperature  along  the  cell  length,  the  density 
will  drop.  In  the  fuel  channel,  the  result  can  be  obtained  based 
on  the  same  analysis. 

Velocity  is  influenced  by  both  the  mass  flow  rate  and  density. 
Owing  to  the  increase  of  density  at  the  entrance,  the  velocity  in 
the  air  channel  will  drop.  With  the  decrease  of  density,  the  mass 
flow  rate  changes  a  little,  so  the  velocity  will  increase  along  the 
cell  length.  While  the  velocity  in  the  fuel  channel  always  rises 
along  the  cell  length,  this  is  mainly  influenced  by  the  increasing 
mass  flow  rate. 

Other  gas  properties,  such  as  specific  heat  capacity,  thermal 
conductivity  and  dynamic  viscosity,  are  primarily  influenced  by 
the  temperature  and  the  gas  composition.  This  has  been  dis¬ 
cussed  by  the  author  Todd  and  Young  [22] .  All  the  profiles  of 
relevant  gas  properties  are  consistent  with  Ref.  [3]. 

As  the  oxygen  ions  transfer  from  the  air  channel  to  the  fuel 
channel,  the  mass  flow  rate  decreases  in  the  air  channel,  while 
increases  in  the  fuel  channel  along  the  cell.  Fig.  8  illustrates  the 
distribution  of  mass  flow  rate  in  both  the  air  and  fuel  channels. 
The  total  mass  in  the  fuel  cell  is  also  provided  to  prove  that  the 
VR  model  accurately  satisfy  the  mass  conservation. 

Fig.  9  illustrates  the  individual  contribution  of  various  poten¬ 
tial  losses,  as  well  as  the  current  density  along  the  cell  length. 
The  potential  losses  and  current  density  distribution  are  related 
to  the  composition  of  the  gas  steam  and  the  temperature  along 
the  cell.  For  high-temperature  SOFC,  the  ohmic  activation  loss 
is  only  of  small  portion,  while  the  cathode  and  anode  activa¬ 
tion  losses  are  the  predominant  source  of  losses.  Along  the  cell 
length,  the  open  circuit  voltage  and  cathode  activation  loss  have 
the  same  trend.  It  increases  near  the  entrance  and  decreases  near 
the  exit,  while  the  anode  activation  loss  is  almost  constant.  Thus, 
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Fig.  9.  Contribution  of  various  potential  losses  and  current  density  along  the  pig  jq  Dynamk  response  of  fuel  utilization, 

cell  length. 


improving  the  operating  temperature  can  effectively  decrease 
the  cathode  activation  loss.  The  current  density  is  a  complicated 
physical  property,  which  is  concerned  with  PEN  temperature, 
partial  pressure  of  fuel  and  fuel  utilization,  etc.  [2].  It  decreases 
at  the  entrance  due  to  the  decreasing  of  PEN  temperature.  With 
the  increasing  of  PEN  temperature,  it  will  rise  in  the  center  of 
fuel  cell.  At  the  exit,  the  current  density  will  decrease  owing  to 
the  lower  partial  pressure  of  fuel. 

The  simulation  result  above  is  consistent  with  Refs. 
[1-3,10,12,17,18,24],  which  verifies  that  the  model  established 
in  this  study  is  useful  and  feasible.  For  the  rated  case,  the  pre¬ 
dicted  values  are  as  following:  the  operating  voltage  is  0.603  V; 
the  average  current  density  is  4802.7  Am-2;  the  power  density 
is  2892.7  Wm-2;  the  fuel  electrical  efficiency  is  46.1%. 

5.3.  Dynamic  simulation  results 

In  this  section,  the  transient  behavior  of  the  fuel  cell  is  inves¬ 
tigated.  Based  on  the  above  basic  steady-state  case,  a  fuel  molar 
rate  step  changes  10%  is  imposed  on  the  fuel  channel  when  the 
fuel  cell  system  has  operated  for  200  s.  Figs.  10-12  presents  the 
dynamic  response  of  the  fuel  cell  performance. 


Fuel  utilization  is  a  very  important  performance  parameter, 
because  it  will  determine  whether  or  not  the  whole  system  is  sta¬ 
bly  operating.  Here,  PID  control  technique  is  adopted  to  keep  the 
fuel  utilization  constant  through  adjusting  the  operating  voltage. 
The  definition  of  the  fuel  utilization  is  presented  in  formula  (25). 
Fig.  10  gives  the  response  of  fuel  utilization.  When  fuel  molar 
flow  rate  is  increased  suddenly,  first  the  consumed  fuel  cannot 
change  at  the  moment,  so  fuel  utilization  has  a  sharp  drop  at  the 
time  200  s.  Thus,  the  operating  voltage  decreases  fast,  and  the 
current  density  will  increase  correspondingly.  The  temperature 
of  the  PEN  structure  will  also  increase  due  to  the  electrochem¬ 
ical  reaction  with  the  increasing  consumed  fuel.  Finally  fuel 
utilization  will  be  back  to  the  constant  0.85  owing  to  the  PID 
controller. 

Fig.  1 1  presents  the  dynamic  response  of  the  performances  of 
the  fuel  cell.  When  the  fuel  molar  rate  is  increased,  the  predicted 
operating  voltage  will  drop  to  0.583  V,  average  current  density 
will  increase  to  5282.9  Am-2,  power  density  will  increase  to 
3078.3  Wm-2,  but  fuel  electrical  efficiency  will  decrease  to 
44.6%. 

Fig.  12  presents  the  dynamic  response  of  the  temperatures 
at  the  exit  of  the  fuel  and  air  channels,  interconnector  and  PEN 


Fig.  11.  (a)  Dynamic  response  of  operating  voltage  of  fuel  cell,  (b)  Dynamic  response  of  average  current  density  of  fuel  cell. 
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Fig.  12.  Temperature  dynamic  response  of  fuel  and  air  channel,  interconnect 
and  PEN  at  the  exit. 


Dimensionless  length  direction 

Fig.  13.  PEN  structure  temperature  profile  comparison  along  the  cell  length. 

structure.  As  can  be  seen,  the  temperature  will  increase  when 
the  fuel  molar  rate  is  increased.  The  maximal  PEN  structure 
temperature  will  increase  to  1284  K. 

According  to  Figs.  11  and  12,  the  inertia  delay  time  of  tem¬ 
perature  and  operating  voltage  is  about  500  s  which  is  mainly 
explainable  in  term  of  the  fuel  cell  great  thermal  capacity.  But 
the  inertia  delay  time  of  current  density  is  only  about  160s.  This 
result  can  be  useful  and  meaningful  for  the  design  of  the  control 
system. 

In  some  sense,  increasing  the  fuel  flow  rate  will  improve  the 
performance  of  the  fuel  cell,  but  some  problems  should  be  con¬ 
sidered  carefully,  such  as  larger  gradient  in  temperature  (Fig.  13) 
and  current  density  (Fig.  14).  Therefore,  the  operation  of  system 
regulation  should  be  appropriate  and  cautious. 

6.  Conclusions 

This  paper  developed  the  mathematical  model  of  DIR-SOFC, 
successfully  introduced  the  volume-resistance  characteristic 
modeling  technique  into  the  SOFC  modeling.  Based  on  this  V-R 
modeling  technique,  the  distributed-lumped  parameter  method 


Dimensionless  length  direction 

Fig.  14.  Current  density  distribution  comparison  along  the  cell  length. 

and  the  modular  modeling  idea,  the  distributed  parameter  SOFC 
simulation  model  was  established.  This  non-iterative  model  can 
satisfy  the  requirement  of  the  quick  dynamic  and  real-time  sim¬ 
ulation. 

This  model  can  predict  the  fuel  cell  distribution  properties, 
such  as  the  temperature,  fuel  compositions,  current  density, 
voltage  losses,  etc.  The  detailed  local  gas  properties  is  also  con¬ 
sidered  along  the  cell  length,  such  as  pressure,  mass  flow  rate, 
velocity,  density,  dynamic  viscosity,  thermal  conductivity  and 
heat  capacity,  etc.  It  is  valuable  and  meaningful  for  the  design 
and  optimization  of  the  fuel  cell  system. 

SOFC  system  has  great  thermal  capacity  indeed,  but  current 
density  delay  time  is  very  short,  only  about  160s.  Dynamic  sim¬ 
ulation  could  provide  reference  for  the  control  system  research 
and  design  in  future. 
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